{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyze the export on hpatches - detector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/yoyee/Documents/deepSfm\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "print(module_path)\n",
    "%matplotlib inline\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# folders\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from imageio import imread\n",
    "def load_as_float(path):\n",
    "    return imread(path).astype(np.float32)/255\n",
    "\n",
    "# folders = ['superpoint_pretrained_hpatches', 'superpoint_coco20_8_170k_hpatches', \\\n",
    "#            'magicpoint_coco_new1_10k_hpatches']\n",
    "# folders = ['superpoint_coco2_3_hpatches-a3']\n",
    "folders = ['superpoint_coco40_16_170k_hpatches_nms4_det0.015']\n",
    "# folders = ['subp32_syn1_0_100k_hpatches']\n",
    "# folders = ['subp32_syn1_0_100k_hpatches']\n",
    "# folders = ['superpoint_pretrained_hpatches', 'superpoint_coco20_5_50k_hpatches']\n",
    "\n",
    "# parameters\n",
    "scale = 10\n",
    "count = 1\n",
    "\n",
    "base_path = '../logs'\n",
    "prediction = 'predictions'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "def getLabels(points_path, img_shape):\n",
    "    pnts = np.load(points_path)['pts']\n",
    "    pnts = pnts.astype(int)\n",
    "    labels = np.zeros((img_shape[0], img_shape[1]))\n",
    "    labels[pnts[:, 1], pnts[:, 0]] = 1\n",
    "    labels_2D = torch.tensor(labels[np.newaxis,:,:], dtype=torch.float32)\n",
    "    return labels_2D, pnts\n",
    "\n",
    "def pnts2img(pnts, img_shape):\n",
    "    pnts = pnts.astype(int)\n",
    "    labels = np.zeros((img_shape[0], img_shape[1]))\n",
    "    labels[pnts[:, 0], pnts[:, 1]] = 1  # pnts (y,x)\n",
    "    labels_2D = torch.tensor(labels[np.newaxis,:,:], dtype=torch.float32)\n",
    "    return labels_2D, pnts\n",
    "\n",
    "def warpLabels(pnts, homography, H, W):\n",
    "    \"\"\"\n",
    "    input: \n",
    "        pnts: numpy\n",
    "        homography: numpy\n",
    "    output:\n",
    "        warped_pnts: numpy\n",
    "    \"\"\"\n",
    "    from utils.utils import warp_points\n",
    "    from utils.utils import filter_points\n",
    "    pnts = torch.tensor(pnts).long()\n",
    "    homography = torch.tensor(homography, dtype=torch.float32)\n",
    "    warped_pnts = warp_points(torch.stack((pnts[:, 0], pnts[:, 1]), dim=1),\n",
    "                                   homography) # check the (x, y)\n",
    "    warped_pnts = filter_points(warped_pnts, torch.tensor([W, H])).round().long()\n",
    "    return warped_pnts.numpy()\n",
    "\n",
    "def img_overlap(img_r, img_g, img_gray):  # img_b repeat\n",
    "    img = np.concatenate((img_gray, img_gray, img_gray), axis=0)\n",
    "    img[0, :, :] += img_r[0, :, :]\n",
    "    img[1, :, :] += img_g[0, :, :]\n",
    "    img[img > 1] = 1\n",
    "    img[img < 0] = 0\n",
    "    img = img.transpose([1,2,0])\n",
    "\n",
    "    return img\n",
    "\n",
    "def printImgPnts(image, pnts):\n",
    "    img_shape = image.shape\n",
    "\n",
    "    print(img_shape)\n",
    "    labels_2D, pnts = pnts2img(keypoints, img_shape)\n",
    "    labels_2D = labels_2D.squeeze().numpy()\n",
    "    print(\"labels \", labels_2D.shape)\n",
    "\n",
    "    img_r = np.zeros_like(image)\n",
    "    overlap = img_overlap(img_r[np.newaxis,:,:], \n",
    "                          labels_2D[np.newaxis,:,:], image[np.newaxis,:,:])\n",
    "    print(overlap.shape)\n",
    "    plt.imshow(overlap)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getNumPoints(data, verbose=False):\n",
    "    # read data\n",
    "    image = data['image']\n",
    "    warped_image = data['warped_image']\n",
    "    real_H = data['homography']\n",
    "    keypoints = data['prob'][:,[1, 0]] # (y, x)\n",
    "    desc = data['desc']\n",
    "    warped_keypoints = data['warped_prob'][:,[1, 0]]\n",
    "    if verbose:\n",
    "        print(\"keypoints shape: \", keypoints.shape)\n",
    "        print(\"desc shape: \", desc.shape)\n",
    "    \n",
    "    # unwarp points, filter\n",
    "    from numpy.linalg import inv\n",
    "    H, W = image.shape\n",
    "    unwarped_keypoints = warpLabels(warped_keypoints, inv(real_H), H, W)\n",
    "    \n",
    "    # return number of points\n",
    "    ## average the points in 2 images\n",
    "    return [keypoints.shape[0], warped_keypoints.shape[0], unwarped_keypoints.shape[0], \n",
    "            (keypoints.shape[0]+unwarped_keypoints.shape[0])/2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### repeatibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_k_best(points, k, thresh=0, verbose=False):\n",
    "    \"\"\" Select the k most probable points (and strip their proba).\n",
    "    points has shape (num_points, 3) where the last coordinate is the proba. \"\"\"\n",
    "    sorted_prob = points\n",
    "    if points.shape[1] > 2:\n",
    "        # add thresh\n",
    "        if verbose: print('num of points: ', points.shape)\n",
    "\n",
    "        sorted_prob = points[points[:,2] > thresh, :]\n",
    "        if verbose: print('num of points after thresh: ', sorted_prob.shape)\n",
    "\n",
    "#             sorted_prob = points[points[:, 2].argsort(), :]\n",
    "        sorted_prob = sorted_prob[sorted_prob[:, 2].argsort(), :]\n",
    "\n",
    "        start = min(k, points.shape[0])\n",
    "        sorted_prob = sorted_prob[-start:, :]\n",
    "    return sorted_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from evaluations.detector_evaluation import compute_repeatability\n",
    "def compute_repeatability(data, keep_k_points=300,\n",
    "                          distance_thresh=3, detection_threshold=0, verbose=False):\n",
    "    \"\"\"\n",
    "    Compute the repeatability. The experiment must contain in its output the prediction\n",
    "    on 2 images, an original image and a warped version of it, plus the homography\n",
    "    linking the 2 images.\n",
    "    \"\"\"\n",
    "    from evaluations.detector_evaluation import warp_keypoints\n",
    "    \n",
    "    def filter_keypoints(points, shape):\n",
    "        \"\"\" Keep only the points whose coordinates are\n",
    "        inside the dimensions of shape. \"\"\"\n",
    "        \"\"\"\n",
    "        points:\n",
    "            numpy (N, (x,y))\n",
    "        shape:\n",
    "            (y, x)\n",
    "        \"\"\"\n",
    "        mask = (points[:, 0] >= 0) & (points[:, 0] < shape[1]) &\\\n",
    "               (points[:, 1] >= 0) & (points[:, 1] < shape[0])\n",
    "        return points[mask, :]\n",
    "\n",
    "    def keep_true_keypoints(points, H, shape):\n",
    "        \"\"\" Keep only the points whose warped coordinates by H\n",
    "        are still inside shape. \"\"\"\n",
    "        \"\"\"\n",
    "        input:\n",
    "            points: numpy (N, (x,y))\n",
    "            shape: (y, x)\n",
    "        return:\n",
    "            points: numpy (N, (x,y))\n",
    "        \"\"\"\n",
    "        # warped_points = warp_keypoints(points[:, [1, 0]], H)\n",
    "        warped_points = warp_keypoints(points[:, [0, 1]], H)\n",
    "        # warped_points[:, [0, 1]] = warped_points[:, [1, 0]]\n",
    "        mask = (warped_points[:, 0] >= 0) & (warped_points[:, 0] < shape[1]) &\\\n",
    "               (warped_points[:, 1] >= 0) & (warped_points[:, 1] < shape[0])\n",
    "        return points[mask, :]\n",
    "\n",
    "    # paths = get_paths(exper_name)\n",
    "    localization_err = -1\n",
    "    repeatability = []\n",
    "    N1s = []\n",
    "    N2s = []\n",
    "    # for path in paths:\n",
    "    # data = np.load(path)\n",
    "    shape = data['image'].shape\n",
    "    H = data['homography']\n",
    "\n",
    "    # Filter out predictions\n",
    "    keypoints = data['prob']\n",
    "    # warped_keypoints = data['warped_prob'][:, :2]\n",
    "    warped_keypoints = data['warped_prob']\n",
    "    warped_keypoints = keep_true_keypoints(warped_keypoints, np.linalg.inv(H),\n",
    "                                           data['image'].shape)\n",
    "\n",
    "    # Warp the original keypoints with the true homography\n",
    "    true_warped_keypoints = keypoints\n",
    "    # true_warped_keypoints[:,:2] = warp_keypoints(keypoints[:, [1, 0]], H)\n",
    "    true_warped_keypoints[:,:2] = warp_keypoints(keypoints[:, :2], H) # make sure the input fits the (x,y)\n",
    "    # true_warped_keypoints = np.stack([true_warped_keypoints[:, 1],\n",
    "    #                                   true_warped_keypoints[:, 0],\n",
    "    #                                   prob], axis=-1)\n",
    "    true_warped_keypoints = filter_keypoints(true_warped_keypoints, shape)\n",
    "\n",
    "    # Keep only the keep_k_points best predictions\n",
    "    warped_keypoints = select_k_best(warped_keypoints, keep_k_points, \n",
    "                                     thresh=detection_threshold, verbose=verbose)\n",
    "    true_warped_keypoints = select_k_best(true_warped_keypoints, keep_k_points, \n",
    "                                     thresh=detection_threshold, verbose=verbose)\n",
    "\n",
    "    # Compute the repeatability\n",
    "    N1 = true_warped_keypoints.shape[0]\n",
    "    N2 = warped_keypoints.shape[0]\n",
    "    N1s.append(N1)\n",
    "    N2s.append(N2)\n",
    "    true_warped_keypoints = np.expand_dims(true_warped_keypoints, 1)\n",
    "    warped_keypoints = np.expand_dims(warped_keypoints, 0)\n",
    "    # shapes are broadcasted to N1 x N2 x 2:\n",
    "    norm = np.linalg.norm(true_warped_keypoints - warped_keypoints,\n",
    "                          ord=None, axis=2)\n",
    "    count1 = 0\n",
    "    count2 = 0\n",
    "    local_err1, local_err2 = None, None\n",
    "    if N2 != 0:\n",
    "        min1 = np.min(norm, axis=1)\n",
    "        count1 = np.sum(min1 <= distance_thresh)\n",
    "        local_err1 = min1[min1 <= distance_thresh]\n",
    "#         print(\"count1: \", count1.shape)\n",
    "#         print(\"local_err1: \", local_err1)        \n",
    "    if N1 != 0:\n",
    "        min2 = np.min(norm, axis=0)\n",
    "        count2 = np.sum(min2 <= distance_thresh)\n",
    "        local_err2 = min2[min2 <= distance_thresh]\n",
    "\n",
    "    if N1 + N2 > 0:\n",
    "        # repeatability.append((count1 + count2) / (N1 + N2))\n",
    "        repeatability = (count1 + count2) / (N1 + N2)\n",
    "    if count1 + count2 > 0:\n",
    "        localization_err = 0\n",
    "        if local_err1 is not None:\n",
    "            localization_err += (local_err1.sum())/ (count1 + count2)\n",
    "        if local_err2 is not None:\n",
    "            localization_err += (local_err2.sum())/ (count1 + count2)\n",
    "    else:\n",
    "        repeatability = 0\n",
    "    if verbose:\n",
    "        print(\"Average number of points in the first image: \" + str(np.mean(N1s)))\n",
    "        print(\"Average number of points in the second image: \" + str(np.mean(N2s)))\n",
    "    # return np.mean(repeatability)\n",
    "    return repeatability, localization_err\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_points(data, thresh=0, verbose=False):\n",
    "    from utils.draw import draw_keypoints\n",
    "    from utils.utils import saveImg\n",
    "    img_2D = data['image']\n",
    "    if verbose: print(\"img_2D max: \", img_2D.max())\n",
    "    for t in thresh:\n",
    "        pts = data['prob']\n",
    "        pts = select_k_best(pts, k=1000, \n",
    "                thresh=t, verbose=verbose)\n",
    "\n",
    "        img_pts = draw_keypoints(img_2D * 255, pts.transpose())\n",
    "        if verbose: print(\"img_pts: \", img_pts.shape)\n",
    "        plt.imshow(img_pts/255)\n",
    "        plt.show()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0,)\n",
      "rep_all:  [0.264305   0.264305   0.264305   0.264305   0.264305   0.264305\n",
      " 0.26504746 0.26540252 0.26665658 0.2681864  0.26970865 0.27227083\n",
      " 0.274579   0.27643696 0.27772441 0.2788896  0.27953152 0.28015486\n",
      " 0.27978154 0.27857416]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEICAYAAADsh6tqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4FeX1wPHvISxhX8OWsArIInvEpa6IFbUCKi4ordaFVmut+1Jaq7S2Kt20tXVp3S0qKEgVROWHolaUQALIHsOWsIUlQAiBLOf3x7zR4Xqz3STMvcn5PE+e3Jl5Z+6ZuTNz7vvOe2dEVTHGGGNiQb2gAzDGGGMqypKWMcaYmGFJyxhjTMywpGWMMSZmWNIyxhgTMyxpGWOMiRl1ImmJSK6I9Aw6DlM1IvKRiNxQ0/PUhGiJIxqJyEYRGRV0HCY8EXlQRF4pY/pKETkrtKyIdHXn3rjqjKdOJC1VbaaqGRUpKyIqIr1qOqZIiUgbEckWkU9Dxp8jImtEJE9EFohItwos63IR+Z+b56Mw0+NE5HcislVEDohIqoi0KmeZJ4jIPBHZJSLf+RGg24n9f0Ui8reqrEdQon1fOVZEpKGIrBaRzBp+n+5um9cPGX+ViGwSkYMiMktE2pSxjCEissTtX0tEZIhv2tlun9snIhvDzLtRRA759t33q3UFj5Hq/gKlqgNU9aMw4ze7c29Rdb5v1Cet0B3U8Ciw2j9CRNoBbwG/BtoAKcDrFVjWHuCvwCOlTH8IOBU4BWgB/BDIL2eZBcAbwPXhJrqduJmqNgM6AoeA6VVcDxOsu4HsIN5YRAYAT+Ptmx2APOAfpZRtCLwNvAK0Bl4E3nbjAQ4Cz+GtT2ku8u3D369C3NVa+6hTVLXMP+BeIAs4AKwFznHjXwB+5yt3FpDpG94I3A+sAvYCzwPxvuk/ANKAHOB/wKCQee8FlgOHgfoVWN6NQDreiXg20Nk3TYFevrifBN516/QFcJybttCVPQjkAleUs23GunXYD3wNjHbjO7sY9riYbvTNEwf80pU/ACwBurhppwKLgX3u/6kh73cq8DnwY+BT3/hJwP98w03xkkHf8j5fV/4G4KOQca3dNjiuIssIs8xe3u5VZplrgAxAKrIewEfADe51J7d/3F3Oe3wE/Bb4zG3v94F2blp393lPArYC24C7fPOOcNs7x037O9CwrH2ljH2i1Djc9JPxjoMcYBlwlm/atW47HQA2AFeXs87XAp8Cf8Q7VjYA54dsk9+598sF/gu0BV51cS8GuruyAvwF2OmmrQBO8C2rB96XqPPxHf/lxPdDYBOwG5iMd2yPctPqAfe5bbcb7wtQGzdts9vmue7vFOD3wH98yz4OOAI0D/O+38c7l4lv3OaSz8g3bhSwMcz838QZwfHwAvBPYI7bZ0YBjdxntBnYATwFNHblzwIy8c4Vu9x7X+1bXlnztgbewfsisde9TnLTHgaK8L585gJ/d+MfB7a4z3gJcLrvvR4EZuB9gTwALAUGh9suruwrIcdX/XDvi3ce/lPIdpoN3F7mtixnQx/vVqSzL4jjfB9CeUnrK6AL3rfmz0rKA0PxDoKT8E7i17jyjXzzprl5G1dgeSPdBzvMfZh/Axb6YglNWrvxTkj18Q7U18KVLWfbjMBLLufiHWiJfHtyXYj3bS8eGOJ2npFu2t14B/7xeCeEwXgnjDZ4O9gPXVwT3HBbN1+c21mG405KvlgeB/4ZEt9XwKUVPKDCJa0z8E6g9wLbgXXAzypxkFYkaf0f8GBF1wOXtPBOlOuASRWI4yO8E2AfoLEbfiTkoJqGlyAHus+q5AAcjpdM6ruyq4HbSttXytknyooj0e2TF7j5znXDCS6u/cDxrmwnYEA563wtXo33Rrff3ISXlMUXSzreCb4l3hfBdXgn0vrAS8Dzrux5eCexVnj7az+gk++93gEuJuT4LyO2/ngnrTPwjtU/A4W+bf4LYBGQ5KY/DUwLPQn6lvc2cG/Ie+QCw8O89+3A3JBx7wB3howrK2ntcPvI+xx94j4NyCljvV9w+8b33Gccj/dlYDbesd8c78vDH/Tb82mh2z6NgDPxkl3JflDWvG2BS4Embtp0YFbIMXFDSHwT3Xz1gTvxjvl4N+1Btz+NBxoAd+F9EWrg2y5lJq1w74t3vGwF6rnhdng15Q5l7UPlNQ8WuQ3WX0QaqOpGVf26nHn8/q6qW1R1D16mneDGTwKeVtUvVLVIVV/Eq1Gd7Jv3CTfvoQos72rgOVVdqqqH8Wpkp4hI91LimqmqX6pqIV7SGlJKubJc797zA1UtVtUsVV0jIl3wdsx7VTVfVdOAfwE/cvPdAPxKVdeqZ5mq7gYuBNar6suqWqiq04A1wEVuvluBL1R1SZhYmuEdEH778HbYSCXhndD64CWJ8cCDInJuFZb5DXet6ky8JpoSFVmP/sAC4Deq+kwF3+55VV3n9qU3+O7n/ZCqHlTVFXg1+AkAqrpEVRe5z2Mj3gn0zDLeJ+w+UYE4JgJzVHWOm+8DvKbRC9z0YuAEEWmsqttUdWUF1nmTqj6r3vWEF/GSXYeQWL5W1X3AXOBrVf3QHRPT8b5Ygneyag70xUt6q1V1G4CIXAzEqerMCsRTYjzwjqoudMfqr936lfgpMFlVM930B4HxZVwmqMy+X9Xj5Gq8E3E3vH1wXsk1XlX9VFXLvN4LvK2qn6lqMd75bhJerWKPqh7AqzVeGTLPr1X1sKp+jNc6dLmISFnzqupuVX1TVfPctIcpe79FVV9x8xWq6p/wzvvH+4osUdUZqlqAl0jjOfp8XWmq+iXe9j/HjboS78vzjrLmKzNpqWo6cBvejrNTRF4Tkc6ViGuL7/UmvGYz8D70O0Ukp+QPrwbVuZR5y1teZzdcEncu3jfVxFLi2u57nYe3M1dWF7xvzqE6AyU7kj/WkljKmm9TyLhNQKLb5rfiNaWEk4t3zcmvBV5VPlIlXxamqOohVV0OvMa3J9Kq+iFebXGDb1xF1uNqvCaeGZV4r/I+77D7lYj0EZF3RGS7iOzHOzG0K+N9Svtsy4ujG3BZyPFwGl6N5iBwBd7JfJuIvCsifct4j++8l6rmuZf+9fafGA6FGW7m5v0/vm3K2Skiz4hICxFpCjyGt19WRmd829ut327f9G7ATN92WI335bkD4VVm36/SceISziGXDP6A1xJxekXmdfz7WQJeTWiJb13fc+NL7HXbp0TJvlnmvCLSRESedp1T9uO1/LQq6zqaiNzlOtPsc8trydH7uv8zK8ZruqxMLijNi3hf2nD/Xy5vhnI7Yqjqf1T1NLydSfE6AoBXVW3iK9oxzOxdfK+74lUFwdsAD6tqK99fE1e7+OatK7G8rS4+ANwB1Rbv5FZTtuA1r4TaCrQREf+3t66+WMqar1vIuJL5RuB9U14lItvxmtFGuJNpHLASr5kR+Gb9j3PjI7Xc/fd/DuE+k0j9iKNrWVCx9XgQryn4P9V4Mbu0/eqfeLXd3qraAu/6gpSxnNI+2/JsAV4OOR6aquojAKo6T1XPxdsH1gDPRvAeEVPVJ1R1OF4ttw9eE3dvvFrHJ26ffAvo5PbJ7mUsbhu+7S0iTfCO1RJb8K6/+bdFvKpmEX7/C91neuLVEtaVUnaQq6mUGETkx4lS9v4QrnyJXXhfDgb41rOleh2USrR2x0CJkn2zvHnvxKslneT22zPc+JJYj9qOInI6cA9wOdDa1Rj3hayb/zOrh9cSs5XKCff5vQKMFZHBeE3Ps8pbSJlJS0SOF5GRItII7wLaIb6tyqcBF4jXBbsjXo0s1M9EJMl1QZ3Mtz3BngV+KiIniaepiFwYcqIPp7TlTQN+7LqzNsL7RvyFa9KprB1ARX7T9W/3nueISD0RSRSRvqq6Be8C9x9EJF5EBuE1G5X8zuFfwG9FpLdb90Ei0hbvAm0f8brv1heRK/BOEu/gNd90x2tOGgI8AKQCQ1zzz0y85qNLRSTeTV8e0jT1HeJ1aY/Ha8eu5+JtAKBeM/AnwGQRaSQi/fCq7++Us0xxy2zohuPdZ+IvcypezXN6yOwVWY8C4DK8az0vuQOoqn7tvp0OwOvkUrJfNce7npTrajc3hcwXuq+E3Scq8P6vABeJyHkln4mInOX29Q4iMtadvA7j1RaKy15c9RGRE91x2gDvi2q+e/+S68sl++QNeNtjCOFbSUrMAH4gIqeJ12tvCkefh54CHhb3UwcRSRCRsW5atntv/zZ/FW/bne620RTgrZCWjhIf4dXabnX79C1u/P+596rn9rsG3qDEuxhLfnP0PfG698eLyN14NZHPyljXUrnayrPAX0SkvXuPRBE5L6ToQ+49T8frvDa9AvM2xztX57hz5W9Clhm63zbHu36WDdQXkQf4bo10uIhcIl4z7W14++KiSq72d86tqpqJ1/HnZeBNPfpyUHha9kXTQcCXeNXnPXgnrJJOGfF4B/d+vG/lt1N678EcvG/VTXzTR7tgS3pmTcf1+CFML50KLO+neE0zJXEm+aaFdsQoqwPJT108OcDl5Wyfi926H8C7sH2eG5/kYtjjYvqpb5444Fd4FzIPuG1Q0rPnNLyL3vvc/9NKed9r8XXE0G8vHq/B21k/wvX+Kif+a9228f+94JueiNfskIvXe+0nFVhm9zDL3BhS5mm8mkW4+UtdD47uPRgPfOg+z3plxPPNPKHbju/2HtwO3OMre4aLJRcvgU/h6A4w39lXytgnSo3DDZ8EfOz2mWy86xdd8WpXH7t9Isctp38FPtfQ/cN/DITG8ruQz30UkO5en+PWJxfvG/6rQLMw73kWFe89eA1er7fSeg/egddT+QDe8fN737xT3PbJAU52465yyzuI1zGjja/8XOCXvuGheMfWIbyOTUND1iF03/3ITRvgtkNJc+Z8INk37+lAbhnr/AK+845vH/493rG1H68p9Fb/9nTbZ5dbvx9WcN7O7jPOxatx/oSjO0Sc4sbvBZ7AOyc955azDa/W5f9MHuTo3oOpwDBfLKFlS+uIcdT7+uaf6MqdXZH9p6Q3UbUT78d5N6jqh9G4PGNcM1ZJL6jCYKMx5lvi3WHiFVVNCjqWmiYiZ+C1NnTTCiSkqP9xsTHGmNrJNTv/AvhXRRIWWNIqk4j8Ur5726FcEZkbdGwVVUr8ua6NPNJlzi1lmb+sztgrEU+1r2O0E5GnSlnnp6IgtqtLia0qHYNMLeOuk+fgNYH/tcLz1VTzoDHGGFPdrKZljDEmZsT0zWjbtWun3bt3DzoMY4yJKUuWLNmlqgnll4w+MZ20unfvTkpKStBhGGNMTBGR0LvvxAxrHjTGGBMzLGkZY4yJGZa0jDHGxAxLWsYYY2KGJS1jjDExw5KWMcaYmGFJyxhjTMyIOGmJyGgRWSsi6SJyX5jpd4jIKhFZLiLzS56P46Y9JiIrxXtS5hMi3kPZRGS4iKxwy/xmvDEmdu3LK+C1Lzfz7MIM3lm+lSWb9rA15xCFRcfssWCmFonox8XiPTH2SeBcvGe+LBaR2aq6ylcsFe95M3kichPeo7mvEO8BgN/De1YXwKfAmXjPf/kncCPwBd5DEUfjPQ/HGBNDCouKWbg+mzeXZPHB6h0cKfxugqonkNC8EZ1aNqZTy3g6toync8vGdGwZT6eW8XRq1Zj2zRvRIM4ahMy3Ir0jxgi8h8RlAIjIa8BYvAc0AqCqC3zlF+E96Au8h32VPNlW8J4SukNEOgEtVHWRW+ZLwDgsaRkTM1Zt3c+bSzN5Oy2LXblHaNO0IVeN6Mr44Ul0adOE7fvy2bbvENv25Xt/OYfYvj+f9TtzWbgum4NHio5angi0b96Iji0b06lFPJ1axbsE15jOLtF1aBFvia0OiTRpJXL0I7Uz8Z68WprrcclHVT8XkQV4T8gU4O+qulpEkt1y/MtMDF2QiEzCe9IsXbt2jTB8Y0x1yT5wmLfTsnhzaRart+2nQZxwTt8OXDIskbOOb0/D+t8mlJaNG3B8x+Zhl6OqHDhcyLYcL7Ft35fP1n35bHdJLj07l0/Wh09sCc0aebWzkJpaJ/faElvtUeP3HhSRiUAyXhMgItIL6If3SHqAD9xzjw5VZHmq+gzwDEBycrI9V8WYAOQXFDF/9U7eWprJR+uyKSpWBie1ZMrYAVw0qDOtmzas9DJFhBbxDWjRsfTEBrA/v8BLaDmVS2ztmjX6pnbmb5Ised2hRfxRCdZEp0iTVhbQxTec5MYdRURGAZOBM1X1sBt9MbBIVXNdmbnAKcDLfJvISl2mMSYYqkrqlhzeXJLJf5dtZX9+IR1bxHPj6T0ZPzyRXu1LTzTVqUV8A1rEN6BPh/ITW0kT5DZfs2RG9kE+S99N7uHCo+YpSWxDu7Ti4qGJnN23PfEN4mp6dUwlRZq0FgO9RaQHXmK5ErjKX0BEhgJPA6NVdadv0mbgRhH5A17z4JnAX1V1m4jsF5GT8Tpi/Aj4W4TxGWOqSVbOIWYuzeStpVlk7DpIfIN6jB7QkUuHJ3Hqce2Iqxd9nXwrktgO5Bd8c21t+75DbM3JJyvnEB+vy+b9VTtoHl+fCwd2YtzQREZ0b0O9KFzPuiiipKWqhSJyCzAPiAOeU9WVIjIFSFHV2cBUoBkw3fVc36yqY4AZwEhgBV6njPdU9b9u0TcDLwCN8a6BWScMYwJw8HAh7321nTeXZvJ5xm5UYUSPNvz0zOM4f2BHmsc3CDrEKmse34DmYRJbUbHyv693MTM1i/8u28pri7fQuWU8Y4YkcvHQxDKbLk3NE9XYvSyUnJys9jwtY6qHqvLlhj28kZLJ3K+2kXekiG5tm3DJ0CQuGZZIlzZNgg7xmMs7UsgHq3YwKzWLhet3UVSs9OvUgouHdmbM4EQ6towPOsSIiMgSVU0OOo5IWNIypo7btu8Qby7JZMaSTDbuzqNZo/r8YFAnLh2eRHK31thv/D27cg/z7vJtzEzNIm1LDiJw6nFtGTskkfNPiK3apyWtgFjSMiYyhwuL+GDVDqanZPLJ+myKFU7p2ZbLkpM4/4RONG5oHRDKsmHXQWalZvF2WhYbd+fRqH49RvXvwMVDEjmjT0LU90K0pBUQS1rGVM5XWfuYnrKFt5dtJSevgM4t4xk/PInxw7vQtW3da/6rKlUlbUsOs1Kz+O/ybew5eIRWTRrwg0GduHhoIsO6RmdN1ZJWQCxpGVO+vQePMCsti+kpmazatp+G9b3ef5clR2/vv1hUUFTMJ+uzmZW6lfdXbSe/oJiubZowbkhnxg5N5LiEZkGH+A1LWgGxpGVMeEXFysL12cxIyeSDVTs4UlTMoKSWXDY8iTGDE2nZJHauv8Si3MOFzPtqO7PSsvgsfRfFCoOSWjJuSCIXDe5MQvNGgcZnSSsglrSMOdqGXQeZsWQLby7JYvv+fNo0bci4IYlclpxEv04tgg6vTtqxP5//LtvKrLQsvsraT1w94bRe7Rg3tDPf79+Rpo1q/MZE32FJKyCWtIzxflM1Z8U2pqdk8uXGPdQTOOv49lyenMTIvh2ivlNAXbJ+xwFmpWUxK3UrWTmHaNwgjvMGdGDc0ERO69WO+sfo/oiWtAJiScvUZXlHCvnjvHW8tngzeUeK6NmuKZcld+GSYYl0aBGbvx+qK4qLlZRNe5mVlsW7y7ex71AB7Zo15KLBnRk3JJFBSS1rtAOHJa2AWNIyddXSzXu54/U0Nu3J45KhSVx1Upeo7almyna4sIiP1mYzKzWL+at3cqSomJ7tmjJuaCLjhiTWSK9OS1oBsaRl6pqComKemL+eJxek06llY/50+WBO7tk26LBMNdl3qIC5K7YxKy2LRRl7ABjerTXjhnTmwkGdaRPB3fPDsaQVEEtapi5Zv+MAt7+RxldZ+xk/PInfXNQ/pu7CYConK+cQs9O2MjM1k3U7cqlfTzjr+ATGDU1kVL8OVboDvSWtgFjSMnVBcbHy/P828uh7a2jWqD6/v3ggo0/oGHRY5hhRVVZvO8DbaVnMSstix/7DNGtUn+tO68Ed5/aJaJmxnLSOfV9LY0yFZeUc4u7py/jf17sZ1a89f7hkUOC/8THHlojQv3ML+nduwT2j+/JFxm5mpmbRsnHdrGVb0jImCqkqs9KyeODtlRQVK49cMpArTuxiHS3quLh6wqm92nFqr3ZBhxIYS1rGRJm9B4/wq1lf8e6KbSR3a82fLh9Mt7ZNgw7LmKhgScuYKLJg7U7unbGcvXlHuGf08fzkjOPs3oDG+ET082sRGS0ia0UkXUTuCzP9DhFZJSLLRWS+iHRz488WkTTfX76IjHPTzhGRpW78pyLSq2qrZkzsyDtSyOSZK/jx84tp3aQhs372PW4+q5clLGNCVLqmJSJxwJPAuUAmsFhEZqvqKl+xVCBZVfNE5CbgMeAKVV0ADHHLaQOkA++7ef4JjFXV1SJyM/Ar4NrIVsuY2OH/ofCkM3pyx7l9qtSd2ZjaLJLmwRFAuqpmAIjIa8BY4Juk5ZJTiUXAxDDLGQ/MVdW8ktmAkjt6tgS2RhCbMTEj9IfC02482X4obEw5IklaicAW33AmcFIZ5a8H5oYZfyXwZ9/wDcAcETkE7AdODrcwEZkETALo2rVrxaM2JorYD4WNiUyN3lJYRCYCycDUkPGdgIHAPN/o24ELVDUJeJ6jE9o3VPUZVU1W1eSEhISaCdyYGlJcrPz70w1c+LdP2ZqTz1MTh/PHywZbwjKmgiKpaWUBXXzDSW7cUURkFDAZOFNVD4dMvhyYqaoFrmwCMFhVv3DTXwfeiyA2Y6LW1pxD3GU/FDamSiJJWouB3iLSAy9ZXQlc5S8gIkOBp4HRqrozzDImAPf7hvcCLUWkj6quw+vksTqC2IyJSgvW7OS219MoKCq2HwobUwWVTlqqWigit+A17cUBz6nqShGZAqSo6my85sBmwHR3YG5W1TEAItIdr6b2ccgybwTeFJFivCR2XVVWzJhoUFhUzF8+XMeTC76mf6cW/OPqYXRvZz8UNiZSdsNcY2rIzgP5/GJaGp9n7GbCiC785qIB1pXdRAW7Ya4x5ihfZOzm59NS2Z9fwB8vG8z44UlBh2RMrWBJy5hqpKo8vTCDqfPW0q1NE166fgR9O7Yof0ZjTIVY0jKmmuzLK+DO6cv4cPUOLhzYiUcuHWhd2Y2pZpa0jKkGKzL3cfN/lrB9Xz6/uag/157a3XoHGlMDLGkZUwWqyn++3MxDs1fRrllDXv/JKQzr2jrosIyptSxpGRMh787sXzEzNYsz+iTw1yuG0KZpw6DDMqZWs6RlTATSd+Zy86tLWL8zlzvO7cMtZ/einj1GxJgaZ0nLmEqavWwr9725nMYN4nj5upM4rXfdffS5MceaJS1jKuhwYREPv7ualz7fRHK31vz9qmF0bBkfdFjG1CmWtIypgMy9efzs1aUsy9zHjaf34J7RfWkQV6MPSTDGhGFJy5hylNzstrhYeWriMEaf0CnokIypsyxpGVMK/81u+3VqwT/tZrfGBM6SljFhZB84zK3TUvk8YzdXJHfhobF2s1tjooElLWNC7DyQz/h/fs7OA/lMHT+Iy5K7lD+TMeaYsKRljM+B/AKufW4x2QcO858bT7a7WxgTZaz7kzHO4cIiJr20hHU7DvDPicMsYRkThSJKWiIyWkTWiki6iNwXZvodIrJKRJaLyHwR6ebGny0iab6/fBEZ56aJiDwsIutEZLWI3Fq1VTOm4oqKldtf9x7YOPWyQZx1fPugQzLGhFHp5kERiQOeBM4FMoHFIjJbVVf5iqUCyaqaJyI3AY8BV6jqAmCIW04bIB14381zLdAF6KuqxSJiZw1zTKgqD/13JXNWbOdXF/bj4qH2wEZjolUkNa0RQLqqZqjqEeA1YKy/gKouUNU8N7gICHcWGA/M9ZW7CZiiqsVuGTsjiM2YSvv7/6Xz0ueb+MkZPbnh9J5Bh2OMKUMkSSsR2OIbznTjSnM9MDfM+CuBab7h44ArRCRFROaKSO9wCxORSa5MSnZ2diVDN+Zo077czJ8+WMclwxK5d3TfoMMxxpSjRjtiiMhEIBmYGjK+EzAQmOcb3QjIV9Vk4FnguXDLVNVnVDVZVZMTEhJqJnBTJ8xbuZ3JM1dw1vEJPHrpILtLuzExIJKklYV37alEkht3FBEZBUwGxqjq4ZDJlwMzVbXANy4TeMu9ngkMiiA2Yyrkyw17uHVaKoOSWvGPq4fZfQSNiRGRHKmLgd4i0kNEGuI18832FxCRocDTeAkr3LWpCRzdNAgwCzjbvT4TWBdBbMaUa832/dzw4mISWzfmuWtPpElD+7miMbGi0kerqhaKyC14TXtxwHOqulJEpgApqjobrzmwGTBdRAA2q+oYABHpjldT+zhk0Y8Ar4rI7UAucENEa2RMGTL35nHNc1/SuGEcL103wp40bEyMEVUNOoaIJScna0pKStBhmBix5+ARxj/1P3YdOMz0n57K8R2bBx2SMYEQkSWu/0DMsXYRUyfkHSnkxy8sJmvvIV6+/iRLWMbEKLv6bGq9gqJibnplKSsyc/jbhKGM6NEm6JCMMRGympap1YqLlXtmLOfjddk8cslAvj+gY9AhGWOqwGpaplZ75L01zEzN4q7v9+HKEV2DDscYU0WWtEyt9czCr3lmYQbXnNKNn53dK+hwjDHVwJKWqZXeWprJ7+es4cKBnXjgogG4n14YY2KcJS1T6yxYu5N7Zizn1OPa8ucrBhNnt2cyptawpGVqldTNe7n5laX06dCcp384nEb144IOyRhTjSxpmVrj6+xcrnthMQnNG/HCdSfSPL5B0CEZY6qZJS1TK2zfl8+P/v0lcfWEl64bQfvm8UGHZIypAZa0TMzbd6iAa577kpy8I7zw4xF0b9c06JCMMTXEflxsYlpRsfLTl5eQsSuXF348ghMSWwYdkjGmBlnSMjHtyQXpfJ6xm8cuHcT3erULOhxjTA2z5kETsxZv3MNfP1zH2CGduSw5KehwjDHHgCUtE5Ny8o7wi2mpJLVuwu/GnWA/HjamjrDmQRNzVJX73lzBzgOHmXHTqda13Zg6JKKaloiMFpG1IpIuIveFmX6HiKwSkeUiMl9EurnxZ4tImu8vX0TGhcz7hIjkRrY6pi549YvNvLdyO3efdzwgGRtqAAAZnElEQVRDurQKOhxjzDFU6aQlInHAk8D5QH9ggoj0DymWCiSr6iBgBvAYgKouUNUhqjoEGAnkAe/7lp0MtI5kRUzdsHb7AX77zipO792OG0/vGXQ4xphjLJKa1gggXVUzVPUI8Bow1l/AJac8N7gICHeVfDwwt6ScS4ZTgXsiiMnUAYeOFPHzaUtpHl+fP10+mHp2T0Fj6pxIklYisMU3nOnGleZ6YG6Y8VcC03zDtwCzVXVbWW8uIpNEJEVEUrKzsysYsqkNfvvuKtbtyOXPlw+xO14YU0fVaEcMEZkIJANnhozvBAwE5rnhzsBlwFnlLVNVnwGeAUhOTtbqjdhEq7krtvGfLzbzkzN6ckafhKDDMcYEJJKklQV08Q0nuXFHEZFRwGTgTFU9HDL5cmCmqha44aFALyDddV1uIiLpqmpP7jNk7s3j3jeXMzipJXd+//igwzHGBCiSpLUY6C0iPfCS1ZXAVf4CIjIUeBoYrao7wyxjAnB/yYCqvgt09M2fawnLABQWFXPba2kUKzwxYSgN69tPC42pyyp9BlDVQrzrT/OA1cAbqrpSRKaIyBhXbCrQDJjuurbPLplfRLrj1dQ+rmLspg54Yv56Ujbt5eGLT6BbW7sRrjF1XUTXtFR1DjAnZNwDvtejyph3I2V33EBVm0USl6ldPv96N39bkM744UmMHVLmLmOMqSOsrcVEpT0Hj3Db66n0aNuUh8YMCDocY0yUsNs4maijqtwzYxl7Dxbw72tOpGkj202NMR6raZmo8+L/NvLh6p3cd35fez6WMeYolrRMVFm5dR+/n7OGkX3b8+PvdQ86HGNMlLGkZaJG3pFCfj4tlVZNGjB1/CB73Igx5jvsYoGJGg/OXsmGXQd59fqTaNusUdDhGGOikNW0TFSYvWwrb6RkcvNZx3Fqr3ZBh2OMiVKWtEzgtuzJY/JbKxjWtRW3jeoTdDjGmChmScsEqqComJ9PSwWBx68cSoM42yWNMaWza1omUH/+YB1pW3J48qphdGnTJOhwjDFRzr7WmsB8un4XT338NRNGdOHCQZ2CDscYEwMsaZlA7Mo9zO1vpHFcQjMe+IHdpskYUzHWPGiOueJi5c43lrHvUAEvXz+Cxg3jgg7JGBMjrKZljrnnPtvAx+uy+fWF/ejbsUXQ4RhjYoglLXNMrcjcx6PvreG8AR2YeHK3oMMxxsQYS1rmmMkvKOIXr6fSrlkjHr3UbtNkjKm8iJKWiIwWkbUiki4i94WZfoeIrBKR5SIyX0S6ufFnuycZl/zli8g4N+1Vt8yvROQ5EWlQtVUz0eYvH64jI/sgU8cPplWThkGHY4yJQZVOWiISBzwJnA/0ByaISP+QYqlAsqoOAmYAjwGo6gJVHaKqQ4CRQB7wvpvnVaAvMBBoDNxQ+dUx0WrZlhyeXZjBlSd24bTedpsmY0xkIqlpjQDSVTVDVY8ArwFj/QVccspzg4uApDDLGQ/MLSmnqnPUAb4sZR4Tgw4XFnH3jGW0bx7PLy/sF3Q4xpgYFknSSgS2+IYz3bjSXA/MDTP+SmBa6EjXLPhD4L1wCxORSSKSIiIp2dnZFQ7aBOfJBV+zbkcuv7/kBFrEW6uvMSZyNdoRQ0QmAsnA1JDxnfCaAeeFme0fwEJV/STcMlX1GVVNVtXkhISE6g7ZVLNVW/fzjwXpXDI0kZF9OwQdjjEmxkXy4+IsoItvOMmNO4qIjAImA2eq6uGQyZcDM1W1IGSe3wAJwE8iiMtEmYKiYu6esYxWTRrywEWhlz2NMabyIqlpLQZ6i0gPEWmI18w3219ARIYCTwNjVHVnmGVMIKRpUERuAM4DJqhqcQRxmSjzzMIMVm7dz+/GDbDegsaYalHppKWqhcAteE17q4E3VHWliEwRkTGu2FSgGTDddW3/JqmJSHe8mtrHIYt+CugAfO7meaCysZnosX7HAR7/cD0XDuzE6BPsZrjGmOoR0b0HVXUOMCdk3AO+16PKmHcjYTpuqKrdB7GWKCpW7p6xnKaN4nhorN0M1xhTfSxRmGr3/GcbSNuSw+NXDqFds0ZBh2OMqUXsNk6mWm3YdZCp89Yyql8HxgzuHHQ4xphaxpKWqTbFxcq9by6nYf16PHzxCXZvQWNMtbOkZarNK19s4ssNe/j1D/rToUV80OEYY2ohS1qmWmzZk8cjc9dweu92XDbc7sBljKkZlrRMlakq97+1AgH+cMlAaxY0xtQYS1qmyt5I2cKn6bu474J+JLVuEnQ4xphazJKWqZLt+/L53TurOblnG64e0TXocIwxtZwlLRMxVWXyzBUUFBfz6KWDqFfPmgWNMTXLkpaJ2Ky0LOav2cnd5/WlW9umQYdjjKkDLGmZiOw8kM+Ds1cxrGsrrj21e9DhGGPqCEtaJiK/eXslhwqKeGz8YOKsWdAYc4xY0jKVNmfFNuZ+tZ3bRvWmV/tmQYdjjKlDLGmZStlz8AgPvP0VAxNbMun0nkGHY4ypY+wu76ZSpvx3JfsOFfDKDSdRP86+8xhjji0765gK+3DVDmalbeVnZ/eib8cWQYdjjKmDIkpaIjJaRNaKSLqI3Bdm+h0iskpElovIfBHp5saf7Z5KXPKXLyLj3LQeIvKFW+brImLPZ48i+w4V8MuZK+jbsTk3n9Ur6HCMMXVUpZOWiMQBTwLnA/2BCSLSP6RYKpCsqoOAGcBjAKq6QFWHqOoQYCSQB7zv5nkU+Iuq9gL2AtdHsD6mhjz87ip2HzzC1PGDaVjfKujGmGBEcvYZAaSraoaqHgFeA8b6C7jklOcGFwHhbvs9Hpirqnni3WF1JF6CA3gRGBdBbKYGLFyXzRspmUw6oycDk1oGHY4xpg6LJGklAlt8w5luXGmuB+aGGX8lMM29bgvkqGphecsUkUkikiIiKdnZ2ZUK3FRe7uFC7n9rBcclNOUX5/QOOhxjTB1Xo+08IjIRSAamhozvBAwE5lV2mar6jKomq2pyQkJC9QRqSvXo3DVs3XeIx8YPJr5BXNDhGGPquEi6vGcBXXzDSW7cUURkFDAZOFNVD4dMvhyYqaoFbng30EpE6rvaVthlmmPr86938/KiTVx/Wg+Gd2sddDjGGBNRTWsx0Nv19muI18w3219ARIYCTwNjVHVnmGVM4NumQVRVgQV417kArgHejiA2U00KioqZPGsFXds04a7vHx90OMYYA0SQtFxN6Ba8pr3VwBuqulJEpojIGFdsKtAMmO66tn+T1ESkO15N7eOQRd8L3CEi6XjXuP5d2dhM9Xnty81kZB/k1z/oT+OG1ixojIkOEd0RQ1XnAHNCxj3gez2qjHk3EqaThapm4PVMNAE7kF/AXz9cz0k92jCqX/ugwzHGmG/YbZzMdzz18dfsPniE5y/sh/drBGOMiQ72K1FzlG37DvGvTzYwdkhnBiW1CjocY4w5iiUtc5Q/zluHKtb5whgTlSxpmW+s3LqPt1Iz+fH3utOlTZOgwzHGmO+wpGUAUFV+P2c1LRs34Oaz7Ya4xpjoZEnLAPDRumw+S9/NrSN707Jxg6DDMcaYsCxpGQqLivnDnNV0a9uEiSd3CzocY4wplSUtw4wlmazbkcu9o/vaY0eMMVHNzlB13MHDhfz5g3UM69qK80/oGHQ4xhhTJktaddyzn2Sw88BhJtsPiY0xMcCSVh22c38+zyzM4IKBHRnerU3Q4RhjTLksadVhf/lwHQVFxdxzXt+gQzHGmAqxpFVHrdtxgNcXb2Hiyd3o3q5p0OEYY0yFWNKqo/4wZzVNG9Xn1pG9gw7FGGMqzJJWHfRZ+i4WrM3mlrN70bppw6DDMcaYCrOkVccUF3u3a0ps1ZhrTu0edDjGGFMpESUtERktImtFJF1E7gsz/Q4RWSUiy0Vkvoh0803rKiLvi8hqV6a7G3+OiCx1Tzr+VETsBng1YFZaFiu37ufu844nvoE9kdgYE1sqnbREJA54Ejgf6A9MEJH+IcVSgWRVHQTMAB7zTXsJmKqq/fCeVLzTjf8ncLWqDgH+A/yqsrGZsuUXFPHHeWsZmNiSMYM7Bx2OMcZUWiQ1rRFAuqpmqOoR4DVgrL+Aqi5Q1Tw3uAhIAnDJrb6qfuDK5frKKdDCvW4JbI0gNlOG5z/byNZ9+dx/QV/q1bMfEhtjYk/9COZJBLb4hjOBk8oofz0w173uA+SIyFtAD+BD4D5VLQJuAOaIyCFgP3ByuIWJyCRgEkDXrl0jCL9u2nPwCP9YkM7Ivu059bh2QYdjjDERqdGOGCIyEUgGprpR9YHTgbuAE4GewLVu2u3ABaqaBDwP/DncMlX1GVVNVtXkhISEGoy+dnli/noOHink/vPth8TGmNgVSdLKArr4hpPcuKOIyChgMjBGVQ+70ZlAmmtaLARmAcNEJAEYrKpfuHKvA6dGEJsJY+Oug7yyaBNXnNiF3h2aBx2OMcZELJKktRjoLSI9RKQhcCUw219ARIYCT+MlrJ0h87ZySQpgJLAK2Au0FJE+bvy5wOoIYjNhTJ23lob163H7qD7lFzbGmChW6WtaqlooIrcA84A44DlVXSkiU4AUVZ2N1xzYDJju7hy+WVXHqGqRiNwFzBdvwhLgWbfMG4E3RaQYL4ldVy1rWMct3byXd1ds4xfn9KZ9i/igwzHGmCoRVQ06hoglJydrSkpK0GFELVXlsqc+Z+PuPD6++yyaNoqk340xprYRkSWqmhx0HJGwO2LUYvNW7iBl015uP7e3JSxjTK1gSauWKigq5tH31tCrfTOuSO5S/gzGGBMDLGnVUtO+3MyGXQe5b3Rf6sfZx2yMqR3sbFYLHcgv4PEP13NSjzac06990OEYY0y1sQsdtdDTH2ew++ARnrugH673pjHG1ApW06pltu/L51+fZnDR4M4M7tIq6HCMMaZaWdKqZf70/lqKi+Ge844POhRjjKl2lrRqkdXb9jNjaSY/OqUbXdo0CTocY4ypdpa0apE/zF1D80b1uWWkPT/TGFM7WdKqJT5Zn83Cddn8fGRvWjVpGHQ4xhhTIyxp1QKHC4v47TurSGzVmB+e0i3ocIwxpsZYl/da4PEP17NuRy7PX3si8Q3igg7HGGNqjNW0Ylzalhye+vhrLk9O4uy+9kNiY0ztZkkrhuUXFHHnG2l0aBHPr37QP+hwjDGmxlnzYAz7y4fr+Dr7IC9eN4IW8Q2CDscYY2qc1bRi1JJNe3l2YQYTRnThzD4J5c9gjDG1QERJS0RGi8haEUkXkfvCTL9DRFaJyHIRmS8i3XzTuorI+yKy2pXp7saLiDwsIuvctFsjXanaLr+giLunL6NTy8b88oJ+QYdjjDHHTKWbB0UkDngSOBfIBBaLyGxVXeUrlgokq2qeiNwEPAZc4aa9BDysqh+ISDOg2I2/FugC9FXVYhGxXgWl+OO8tWTsOsgr159Ec2sWNMbUIZHUtEYA6aqaoapHgNeAsf4CqrpAVfPc4CIgCUBE+gP1VfUDVy7XV+4mYIqqFrtpOyOIrdZbvHEP//5sAxNP7sppvdsFHY4xxhxTkSStRGCLbzjTjSvN9cBc97oPkCMib4lIqohMdTU3gOOAK0QkRUTmikjvcAsTkUmuTEp2dnYE4ceuQ0e8ZsHEVo25/3xrFjTG1D012hFDRCYCycBUN6o+cDpwF3Ai0BOvWRCgEZCvqsnAs8Bz4Zapqs+oarKqJick1K0OCI/NW8PG3XlMHT+Ypo2s46cxpu6JJGll4V17KpHkxh1FREYBk4ExqnrYjc4E0lzTYiEwCxjmm/aWez0TGBRBbLXWoozdPP/ZRq45pRunHNc26HCMMSYQkSStxUBvEekhIg2BK4HZ/gIiMhR4Gi9h7QyZt5WIlFSRRgIlHThmAWe712cC6yKIrVY6eLiQe2Ysp2ubJtx7ft+gwzHGmMBUuo1JVQtF5BZgHhAHPKeqK0VkCpCiqrPxmgObAdPd4943q+oYVS0SkbuA+eJNWILXFAjwCPCqiNwO5AI3VHXlaotH31vDlr15vD7pFJo0tGZBY0zdFdEZUFXnAHNCxj3gez2qjHk/IEzTn6rmABdGEk9t9r/0Xbz0+Sau+14PRvRoE3Q4xhgTKLsjRhTbn1/APW8up0e7ptx93vFBh2OMMYGztqYodehIETe8kML2ffm8/pOTadzQHjlijDFW04pChwuLmPRyCimb9vCXK4YwvJs1CxpjDFhNK+oUFhVz67RUPlm/i8cuHcRFgzsHHZIxxkQNq2lFkeJi5Z4Zy5m3cgcP/KA/l5/YpfyZjDGmDrGkFUUen7+et1KzuPPcPlx3Wo+gwzHGmKhjSStKLM/M4e8L0rl4aCK3jOwVdDjGGBOVLGlFgfyCIu54YxkJzRrx4EUDcD/INsYYE8I6YkSBP3+wjvSdubx43QhaNrHnYxljTGmsphWwxRv38OwnGVx1UlfO7FO37lpvjDGVZUkrQAcPF3LnG8tIat2YX15gz8cyxpjyWPNggB6Z690Id9qNJ9PMno9ljDHlsppWQD5Zn83Li7wb4Z7c056PZYwxFWFJKwD78wu4Z8ZyeibYjXCNMaYy6mSb1FdZ+7j99bTA3v9AfiE7D+Tz1s3fI76B3QjXGGMqqk4mrfgGcfTu0CzQGL7fvyNDurQKNAZjjIk1ESUtERkNPI735OJ/qeojIdPvwHvycCGQDVynqpvctK7Av4AugAIXqOpG37xPuPI1llV6tW/GP64eXlOLN8YYU0MqfU1LROKAJ4Hzgf7ABBHpH1IsFUhW1UHADOAx37SXgKmq2g8YAez0LTsZaF3ZmIwxxtQNkXTEGAGkq2qGqh4BXgPG+guo6gJVzXODi4AkAJfc6qvqB65cbkk5lwynAvdEtCbGGGNqvUiSViKwxTec6caV5npgrnvdB8gRkbdEJFVEprpkBXALMFtVt5X15iIySURSRCQlOzs7gvCNMcbEqhrtiCEiE4Fk4Ezf+50ODAU2A68D14rIXOAy4KzylqmqzwDPACQnJ2v1R22MMSZaRZK0svA6UZRIcuOOIiKjgMnAmap62I3OBNJUNcOVmQWcDGwHegHp7g7nTUQkXVXtGR3GGGO+EUnSWgz0FpEeeMnqSuAqfwERGQo8DYxW1Z0h87YSkQRVzQZGAimq+i7Q0Td/riUsY4wxoSp9TUtVC/GuP80DVgNvqOpKEZkiImNcsalAM2C6iKSJyGw3bxFwFzBfRFYAAjxbDethjDGmDhDV2L0slJycrCkpKUGHYYwxMUVElqhqctBxRCKmk5aIZAObqrCIdsCuagqnpsRCjGBxVjeLs3pZnEfrpqox+QC/mE5aVSUiKdH+bSMWYgSLs7pZnNXL4qw97C7vxhhjYoYlLWOMMTGjrietZ4IOoAJiIUawOKubxVm9LM5aok5f0zLGGBNb6npNyxhjTAyxpGWMMSZm1MqkJSKjRWStiKSLyH1hpjcSkdfd9C9EpLtv2v1u/FoROS8a4xSRc0VkiYiscP9HRmOcvuldRSRXRO6K1jhFZJCIfC4iK912jY+2OEWkgYi86OJbLSL311SMFYzzDBFZKiKFIjI+ZNo1IrLe/V0TjXGKyBDfZ75cRK6Ixjh901uISKaI/L0m44x6qlqr/vCepvw10BNoCCwD+oeUuRl4yr2+Enjdve7vyjcCerjlxEVhnEOBzu71CUBWNG5P3/QZwHTgrmiME+8enMuBwW64bZR+7lcBr7nXTYCNQPcA4+wODMJ7sOt43/g2QIb739q9bh2FcfYBervXnYFtQKtoi9M3/XHgP8DfayLGWPmrjTWtch9S6YZfdK9nAOeId3v5sXgnhcOqugFId8uLqjhVNVVVt7rxK4HGItIo2uIEEJFxwAYXZ02qSpzfB5ar6jIAVd2t3n0yoy1OBZqKSH2gMXAE2B9UnKq6UVWXA8Uh854HfKCqe1R1L/ABMDra4lTVdaq63r3eivcU9Zq6S0RVticiMhzoALxfQ/HFjNqYtCrykMpvyqh3A+B9eN+uK/uAy6Di9LsUWKrfPv4lauIUkWbAvcBDNRRbtcSJ941bRWSea56pyadnVyXOGcBBvBrBZuCPqronwDhrYt7Kqpb3EpEReDWgr6sprlARxyki9YA/4d1svM6r0YdAmpolIgOAR/FqCtHoQeAvqprrKl7Rqj5wGnAikIf3FIIlqjo/2LC+YwRQhNeU1Rr4REQ+VPd8OhMZEekEvAxco6rfqeVEgZuBOaqaGeXH0TFRG2taFXlI5TdlXFNLS2B3BeeNhjgRkSRgJvAjVa2pb4dVjfMk4DER2QjcBvxSRG6JwjgzgYWquktV84A5wLAojPMq4D1VLVDvOXWf4T0ZPKg4a2LeyqrSe4lIC+BdYLKqLqrm2PyqEucpwC3uOPoj8CMReaR6w4shQV9Uq+4/vG/NGXgdKUoueA4IKfMzjr7Q/YZ7PYCjO2JkUHMX5KsSZytX/pJo3p4hZR6kZjtiVGV7tgaW4nVuqA98CFwYhXHeCzzvXjcFVgGDgorTV/YFvtsRY4Pbrq3d6zZRGGdDYD5wW03tl9URZ8i0a6njHTECD6CGdpALgHV47dOT3bgpwBj3Oh6vN1s68CXQ0zfvZDffWuD8aIwT+BXetY0031/7aIszZBkPUoNJqxo+94l4nUW+Ah6LxjhxD1Z1ca4C7g44zhPxaqkH8WqCK33zXufiTwd+HI1xus+8IOQ4GhJtcYYs41rqeNKy2zgZY4yJGbXxmpYxxphaypKWMcaYmGFJyxhjTMywpGWMMSZmWNIyxhgTMyxpGWOMiRmWtIwxxsSM/wdxQjZv5Ud6uAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# images_num = 10\n",
    "# rows = images_num\n",
    "# cols = len(folders)\n",
    "# plt.figure(figsize=(cols*2*scale, rows*scale))\n",
    "# task_folder = ['repeatibility3']\n",
    "count = 1\n",
    "max_iter = 600\n",
    "verbose = False\n",
    "rep_thd = 1\n",
    "rep_thd = np.array(rep_thd)[...,np.newaxis]\n",
    "det_thd = [0, 0.003, 0.006, 0.009, 0.012, 0.015, 0.018, 0.02, 0.025, 0.030, 0.035, \n",
    "           0.045, 0.055, 0.065, 0.075, 0.085, 0.095, 0.105, 0.125, 0.145]\n",
    "# det_thd = [0.012, 0.015]\n",
    "rep_all = []\n",
    "plot_images = False\n",
    "if verbose: print(\"rep_thd: \", rep_thd)\n",
    "if verbose: print(\"det_thd: \", det_thd)\n",
    "\n",
    "all_num_points = []\n",
    "for folder in folders:\n",
    "    num_points = [] # [[keypoints, warped_keypoints, unwarped_keypoints]]\n",
    "    exp_path = Path(base_path, folder, prediction)\n",
    "    files = os.listdir( exp_path )\n",
    "#     for f in files:\n",
    "    for f in files:\n",
    "    #     print(\"file: \", f)\n",
    "        if count > max_iter:\n",
    "            break\n",
    "        if f[-3:] == 'npz':\n",
    "            if verbose: print(\"file: \", f)\n",
    "            data = np.load(exp_path/f)\n",
    "            if verbose: print(list(data))\n",
    "            reps = []\n",
    "            for t in det_thd:\n",
    "                if verbose: print(\"det thd: \", t)\n",
    "                rep, local_err = compute_repeatability(data, keep_k_points=300, \n",
    "                                                       distance_thresh=rep_thd, detection_threshold=t, \n",
    "                                                       verbose=verbose)\n",
    "                reps.append(rep)\n",
    "            rep_all.append(reps)\n",
    "            \n",
    "            # plot images\n",
    "            if plot_images:\n",
    "                visualize_points(data, thresh=det_thd, verbose=True)\n",
    "            \n",
    "#             print(\"reps: \", reps)\n",
    "#             num_points.append(getNumPoints(data))\n",
    "    #         print(\"points: \", getNumPoints(data))\n",
    "            count += 1\n",
    "    # file =  exp_path / '0.npz'\n",
    "\n",
    "    num_points = np.array(num_points)\n",
    "    print(num_points.shape)\n",
    "    all_num_points.append(num_points)\n",
    "\n",
    "    rep_all = np.array(rep_all)\n",
    "    print('rep_all: ', rep_all.mean(axis=0))\n",
    "    \n",
    "    plt.plot(det_thd, rep_all.mean(axis=0))\n",
    "    plt.title(folder + \": repeatability\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py36_pytorch",
   "language": "python",
   "name": "py36_pytorch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
